Master File¶
This notebook runs the GP backtest, NN-IV training and backtesting, SSVI backtesting and trinomial tree backtesting. Finally all results are compared.
Help on Importing files¶
If you run this notebook on google colab you need to upload python scripts on the left panel. To that end click on the left "Files" (or "Fichiers" in French) and drag and drop :
- python scripts from the "code" folder of github repository.
- csv files or xls files from "data" folder.
- Tensorflow model (files with .data or .index extensions) in the "Results" folder if you want to use the neural network while avoiding the training step.
- csv files in "Results" folder if you want to load Results presented in our paper.
In [2]:
import sys
formerPath = sys.path
sys.path.append('./code/')
sys.path.append('./BS/')
In [24]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sklearn as skl
import bootstrapping
import dataSetConstruction
import backtest
import BS
import loadData
import plotTools
import SSVI
import importlib
import neuralNetwork
from Bisect import bisect
from newton import newton
from BSImplVol import bsimpvol
%matplotlib inline
Load data¶
In order to reproduce our paper experiments, execute cells from part "Load preformatted data".
Each source of data produces the following objects :
- bootstrap manages discounting and dividend.
- dataSet contains the training set.
- dataSetTest contains the testing set.
- $S0$ the spot value of the underlying
In [4]:
importlib.reload(loadData)
Out[4]:
In [5]:
pd.options.mode.chained_assignment = 'raise'
In [6]:
workingFolder = "./data/"
fileName = "Option_SPX_18_Mai_2019Feuille2.xlsm"
asOfDate = "2019-05-18"
trainingSet, testingSet, bootstrap, S0 = loadData.loadCBOTData(workingFolder, fileName, asOfDate)
In [6]:
importlib.reload(bootstrapping)
Out[6]:
In [7]:
workingFolder = "./data/"
trainingSet, testingSet, bootstrap, S0 = loadData.loadFormattedData(workingFolder)
Formatting data¶
Boostsrapping Rate Curve¶
We assume a piecewise constant discount short rate $r$ and a piecewise constant dividend short rate $q$.
We estimate the "zero coupon dividend" $D(T) = e^{-\int_{0}^{T} q_s ds}$ by regressing it against maturity : $$e^{-\int_{0}^{T} q_s ds} = \frac{C(T,K) - P(T,K) + K e^{-\int_{0}^{T} r_s ds}}{S_0}$$
Then we have $\hat{q}_t = - \frac{ \log{D(\overline{T})} - \log{D(\underline{T})} }{ \overline{T} - \underline{T} }$ with $\overline{T}$ the smallest discretized maturity greater than $T$ and $\underline{T}$ the grestest discretized maturity inferior than $T$.
bootstrap object has several members :
- riskFreeIntegral corresponds to $I_T = \int_{0}^{T} r_u du$.
- riskFreeSpline corresponds to $r_u$ evaluated on a subgrid. Interpolated as step function from zero coupons
- divSpreadIntegral corresponds to $I_T = \int_{0}^{T} q_u du$, can be negative.
- divSpline corresponds to $q_u$ evaluated on a subgrid, can be negative.
These curve should satisfy the call-put parity.
Change of variable¶
Neural network on modified prices with modified strike as input such that discounting and dividend don't intervene in Dupire formula calculation.
- In presence of dividend rate $d$ and risk free rate $r$ Dupire formula is : $$\sigma^2(T,K) = 2 \frac{ \partial_T P(T,K) + (r-q) K \partial_K P(T,K) + qP(T,K)}{K² \partial_{K}^2 P(T,K)}$$ with Strike $K$, Maturity $T$, dividend rate $q$ and risk-free rate $r$, $P$ our pricing function.
- We apply the following change of variable : $$ w(T,k) = \exp{(\int_{0}^{T} q_t dt)} P(T,K)$$ with $K = k \exp{(\int_{0}^{T} (r_t - q_t) dt)} $.
- Then Dupire equation becomes : $\sigma^2(T,K) = 2 \frac{ \partial_T w(T,k)}{k² \partial_{k}^2 w(T,k)}$.
- If we learn the mapping $v$ with a neural network then we obtain quickly by adjoint differentiation $\partial_T w$ and $\partial_{k²}^2 w$ and therefore $\sigma$.
- $k$ corresponds to "ChangedStrike" column.
- $\exp{(\int_{0}^{T} q_t dt)}$ corresponds to "DividendFactor" column.
In [8]:
dataSet = trainingSet #Training set
dataSetTest = testingSet #Testing set
dataSet.head()
Out[8]:
In [9]:
dataSetTest.head()
Out[9]:
In [10]:
dfCurve = dataSetConstruction.savingData(bootstrap,
dataSet,
dataSetTest,
workingFolder)
In [11]:
threshold = pd.concat([dataSet, dataSetTest])["Price"].sort_values().iloc[int(0.05 * pd.concat([dataSet, dataSetTest]).shape[0])]
threshold
Out[11]:
In [15]:
dataSet.shape
Out[15]:
In [16]:
dataSetTest.shape
Out[16]:
In [17]:
S0
Out[17]:
In [12]:
def getCall(df):
return df[df["OptionType"]==1]
def getPut(df):
return df[df["OptionType"]==-1]
In [19]:
plotTools.plotSerie(getPut(trainingSet)["ImpliedVol"],
yMin=0,
yMax=3*S0,
Title = 'True Price Surface')
In [20]:
plotTools.plotSerie(getPut(testingSet)["ImpliedVol"],
yMin=0,
yMax=3*S0,
Title = 'True Price Surface')
In [21]:
plotTools.plot2Series(getPut(testingSet)["Price"],
getPut(trainingSet)["Price"],
yMin=0,
yMax=3*S0,
Title = 'True Price Surface')
Scaling methods¶
Use min-max of scaling strike between 0 et 1 for improving stability of neural network training.
In [13]:
activateScaling = False
transformCustom = dataSetConstruction.transformCustomMinMax if activateScaling else dataSetConstruction.transformCustomId
inverseTransform = dataSetConstruction.inverseTransformMinMax if activateScaling else dataSetConstruction.inverseTransformId
inverseTransformColumn = dataSetConstruction.inverseTransformColumnMinMax if activateScaling else dataSetConstruction.inverseTransformColumnId
inverseTransformColumnGreeks = dataSetConstruction.inverseTransformColumnGreeksMinMax if activateScaling else dataSetConstruction.inverseTransformColumnGreeksId
In [14]:
scaler = skl.preprocessing.MinMaxScaler(feature_range=(0, 1))
scaler.fit(dataSet)
scaledDataSet = transformCustom(dataSet, scaler)
scaledDataSetTest = transformCustom(dataSetTest, scaler)
In [24]:
scaledDataSet.head()
Out[24]:
In [15]:
#Search strike for ATM option
midS0 = dataSet[dataSet.index.get_level_values("Strike") >= S0].index.get_level_values("Strike").min()
In [17]:
KMin = 0
KMax = 2.0 * S0
GP Backtesting¶
The purpose of this section is to load the GP local volatility surface and perform the Monte-Carlo backtest of the option prices uses the GP local volatility surface. Note that the GP local volatility is generated by running the Matlab code in the "code/GP" folder. See Section "GP Local Volatility Backtests" below for further details of the backtests.
Load GP Results¶
This section loads result from the Matlab experiments. See code/GP folder to access the matlab script.
In [28]:
pathFolder = "./data/"
fileName = "GP_output_Put_Price_testing_set.xlsx"
putTestingGP = pd.read_excel(pathFolder + fileName,
header=0,
sheet_name = "Sheet1")
putTestingGP["Strike"] = np.exp(+ bootstrap.discountIntegral(putTestingGP["T"])
- bootstrap.dividendIntegral(putTestingGP["T"])) * putTestingGP["K"]
putTestingGP["Maturity"] = putTestingGP["T"]
putTestingGP = putTestingGP.set_index(["Strike","Maturity"], drop=False)
putTestingGP = pd.DataFrame(putTestingGP.values, columns = putTestingGP.columns, index = dataSetTest.index)
putTestingGP["Strike"] = dataSetTest["Strike"]
putTestingGP["Maturity"] = dataSetTest["Maturity"]
putTestingGP.head()
Out[28]:
In [29]:
putTestingGP.tail()
Out[29]:
In [30]:
dataSet.head()
Out[30]:
In [31]:
ImpVolPutTesting = BS.vectorizedImpliedVolatilityCalibration(S0, bootstrap,
putTestingGP["T"],
putTestingGP["Strike"],
-1 * np.ones_like(putTestingGP["Strike"]),
putTestingGP["GP_Put_price"],
removeNaN= False)
ImpVolPutTesting = pd.Series(ImpVolPutTesting, index = putTestingGP.set_index(["Strike","Maturity"], drop=False).index).sort_index()
In [32]:
plotTools.predictionDiagnosis(putTestingGP["GP_Put_price"],
dataSetTest['Price'],#ImpVolPutTestingRef, #
" Price ",
yMin=2400,
yMax=3600,
az = 225,
threshold=None)
In [33]:
keptPrices = ImpVolPutTesting.dropna().index
plotTools.predictionDiagnosis(plotTools.selectIndex(ImpVolPutTesting, keptPrices),
plotTools.selectIndex(dataSetTest["ImpliedVol"], keptPrices),
" Implied vol ",
yMin=2400,
yMax=3600,
az = 225)
In [34]:
pathFolder = "./data/"
fileName = "GP_output_Put_Price_training_set.xlsx"
putTrainingGP = pd.read_excel(pathFolder + fileName,
header=0,
sheet_name = "Sheet1")
putTrainingGP["Strike"] = np.exp(+ bootstrap.discountIntegral(putTrainingGP["T"])
- bootstrap.dividendIntegral(putTrainingGP["T"])) * putTrainingGP["K"]
putTrainingGP["Maturity"] = putTrainingGP["T"]
putTrainingGP = putTrainingGP.set_index(["Strike","Maturity"], drop=False)
putTrainingGP = pd.DataFrame(putTrainingGP.values, columns = putTrainingGP.columns, index = dataSet.index)
putTrainingGP["StrikeGap"] = dataSet["Strike"] - putTrainingGP["Strike"]
putTrainingGP["Strike"] = dataSet["Strike"]
putTrainingGP["Maturity"] = dataSet["Maturity"]
putTrainingGP["OriginalPrice"] = putTrainingGP["GP_Put_price"] * np.exp(- bootstrap.dividendIntegral(putTrainingGP["T"]))
putTrainingGP.head()
Out[34]:
In [35]:
ImpVolPutTraining = BS.vectorizedImpliedVolatilityCalibration(S0, bootstrap,
putTrainingGP["T"],
putTrainingGP["Strike"],
-1 * np.ones_like(putTrainingGP["Strike"]),
putTrainingGP["GP_Put_price"],
removeNaN= False)
ImpVolPutTraining = pd.Series(ImpVolPutTraining, index = putTrainingGP.set_index(["Strike","Maturity"], drop=False).index).sort_index()
In [36]:
plotTools.predictionDiagnosis(putTrainingGP["GP_Put_price"],
dataSet['Price'],
" Implied vol ",
yMin=2400,
yMax=3600,
az = 30,
threshold=None)
In [37]:
keptPrices = ImpVolPutTraining.dropna().index
plotTools.predictionDiagnosis(plotTools.selectIndex(ImpVolPutTraining, keptPrices),
plotTools.selectIndex(dataSet["ImpliedVol"], keptPrices),
" Implied vol ",
yMin=2400,
yMax=3600,
az = 180)
Load GP local volatility¶
In [38]:
strikeLow = min(dataSet["Strike"].min(),dataSetTest["Strike"].min()) #dataSet["Strike"].min()
strikeUp = max(dataSet["Strike"].max(),dataSetTest["Strike"].max()) #dataSet["Strike"].max()
strikeGrid = np.linspace(strikeLow, strikeUp, 100)
matLow = min(dataSet["Maturity"].min(),dataSetTest["Maturity"].min()) #dataSet["Maturity"].min()
matUp = max(dataSet["Maturity"].max(),dataSetTest["Maturity"].max()) #dataSetTest["Maturity"].max()
matGrid = np.linspace(matLow, matUp, 100)
volLocaleGrid = np.meshgrid(strikeGrid, matGrid)
volLocaleGridDf = pd.DataFrame(np.vstack((np.ravel(volLocaleGrid[0]), np.ravel(volLocaleGrid[1]))).T,
columns = ["Strike", "Maturity"]).set_index(["Strike","Maturity"],drop=False)
volLocaleGridDf["ChangedStrike"] = bootstrap.changeOfVariable(volLocaleGridDf["Strike"], volLocaleGridDf["Maturity"])[0]
volLocaleGridDf["logMoneyness"] = np.log(volLocaleGridDf["ChangedStrike"] / S0)
volLocaleGridDf["OptionType"] = np.ones_like(volLocaleGridDf["logMoneyness"])
In [39]:
workingFolder = "./data/volLocaleAresky/"
#filename = "local_vol_nx_10_nt_27.xlsx" RMSE : 10.073967351737087
#filename = "local_vol_nx_12_nt_27.xlsx" RMSE : 7.758802111118254
#filename = "local_vol_nx_15_nt_27.xlsx" RMSE : 6.252799135416868
#filename = "local_vol_nx_18_nt_27.xlsx" RMSE : 5.083806059940602
#filename = "local_vol_nx_20_nt_27.xlsx" RMSE : 5.050554384554841
#filename = "local_vol_nx_25_nt_27.xlsx" RMSE : 5.499835852015688
filename = "local_vol_nx_20_nt_27.xlsx"
In [40]:
locVolGP = loadData.loadGPLocVol(workingFolder, filename, bootstrap, S0)
In [41]:
locVolGP.head()
Out[41]:
In [42]:
filterLocVol = locVolGP[(locVolGP["Maturity"] <= matUp) & (locVolGP["Maturity"] >= matLow)]
plotTools.plotSerie(filterLocVol["LocalVolatility"],
Title = 'GP local volatility',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [43]:
filterLocVol = locVolGP[(locVolGP["Maturity"] <= matUp) & (locVolGP["Maturity"] >= matLow)]
plotTools.plotSerie(plotTools.convertToLogMoneyness(filterLocVol["LocalVolatility"], S0),
Title = 'GP local volatility',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [44]:
nnGP = lambda x,y : backtest.interpolatedMCLocalVolatility(locVolGP["LocalVolatility"], x, y)
filterLocVol = locVolGP[(locVolGP["Maturity"] <= matUp) & (locVolGP["Maturity"] >= matLow)]
plotTools.plotSerie(nnGP(dataSetTest["Strike"], dataSetTest["Maturity"]),
Title = 'GP local volatility',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [45]:
filterLocVol = locVolGP[(locVolGP["Maturity"] <= matUp) & (locVolGP["Maturity"] >= matLow)]
plotTools.plotSerie(plotTools.convertToLogMoneyness(nnGP(dataSetTest["Strike"], dataSetTest["Maturity"]), S0),
Title = 'GP local volatility',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [46]:
logMin = np.log(KMin/S0),
logMax = 0.1 #np.log(KMax/S0),
plotTools.plot2Series(plotTools.convertToLogMoneyness(nnGP(dataSet["Strike"], dataSet["Maturity"]), S0),
plotTools.convertToLogMoneyness(nnGP(dataSetTest["Strike"], dataSetTest["Maturity"]), S0),
yMin=logMin,
yMax=logMax,
az = 340,
Title = 'True Implied Vol Surface')
In [47]:
filterLocVol = locVolGP[(locVolGP["Maturity"] <= matUp) & (locVolGP["Maturity"] >= matLow)]
plotTools.plotSerie(nnGP(volLocaleGridDf["Strike"], volLocaleGridDf["Maturity"]),
Title = 'GP local volatility',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [77]:
plotTools.saveDataModel(plotTools.removeDuplicateIndex(putTrainingGP["GP_Put_price"]),
plotTools.removeDuplicateIndex(nnGP(dataSet["Strike"], dataSet["Maturity"])),
plotTools.removeDuplicateIndex(ImpVolPutTraining) ,
"./Results/GPTraining")
In [78]:
plotTools.saveDataModel(plotTools.removeDuplicateIndex(putTestingGP["GP_Put_price"]),
plotTools.removeDuplicateIndex(nnGP(dataSetTest["Strike"], dataSetTest["Maturity"])),
plotTools.removeDuplicateIndex(ImpVolPutTesting) ,
"./Results/GPTesting")
GP Local Volatility Backtests¶
During Monte Carlo backtest, each option in testing set is priced with an underlying which is diffused with the following SDE : $$ dS_t = \left( r_t - q_t - \frac{\sigma_{NN}^2(t, S_t)}{2} \right) dt + \sigma_{NN}(t, S_t) dW_t$$ with $\sigma_{NN}$ the neural local volatility function.
Due to computation time issue we avoid to make millions of call to neural network and we interpolate linearly neural local volatility obtained on one the two possible grid :
- the testing grid i.e. nodes $(T,K)$ of the testing set.
- an artificial grid of 10000 points to check local volatility is correctly interpolated/extrapolated. That grid is the smallest rectangle containing the minimum and maximum maturities and the minimum and maximum strikes of our dataset (union of testing and training set).
During PDE backtest, we used a crank-nicholson scheme to revaluate each option in our testing set. Time step corresponds to one day and space grid has 100 points.
In [75]:
nbTimeStep = 100
nbPaths = 100000
In [76]:
importlib.reload(backtest)
Out[76]:
In [77]:
nnGP = lambda x,y : backtest.interpolatedMCLocalVolatility(locVolGP["LocalVolatility"], x, y)
In [78]:
strikeLow = min(dataSet["Strike"].min(),dataSetTest["Strike"].min()) #dataSet["Strike"].min()
strikeUp = max(dataSet["Strike"].max(),dataSetTest["Strike"].max()) #dataSet["Strike"].max()
strikeGrid = np.linspace(strikeLow, strikeUp, 100)
matLow = min(dataSet["Maturity"].min(),dataSetTest["Maturity"].min()) #dataSet["Maturity"].min()
matUp = max(dataSet["Maturity"].max(),dataSetTest["Maturity"].max()) #dataSetTest["Maturity"].max()
matGrid = np.linspace(matLow, matUp, 100)
volLocaleGrid = np.meshgrid(strikeGrid, matGrid)
volLocaleGridDf = pd.DataFrame(np.vstack((np.ravel(volLocaleGrid[0]), np.ravel(volLocaleGrid[1]))).T,
columns = ["Strike", "Maturity"]).set_index(["Strike","Maturity"],drop=False)
volLocaleGridDf["ChangedStrike"] = bootstrap.changeOfVariable(volLocaleGridDf["Strike"], volLocaleGridDf["Maturity"])[0]
volLocaleGridDf["logMoneyness"] = np.log(volLocaleGridDf["ChangedStrike"] / S0)
volLocaleGridDf["OptionType"] = np.ones_like(volLocaleGridDf["logMoneyness"])
In [79]:
volLocaleGridDf.head()
Out[79]:
In [86]:
mcResGPTest = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnGP)
mcResGPTest.head()
Out[86]:
In [87]:
mcResGPTest.to_csv(workingFolder + "mcResGPTest.csv")
In [88]:
plotTools.predictionDiagnosis(mcResGPTest["Price"],
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [89]:
plotTools.plotSerie(mcResGPTest["stdPrice"],
Title = 'GP backtested price std',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [90]:
testingDataSetGP = dataSetConstruction.interpolatedLocalVolatility(locVolGP,
dataSetTest["Price"])
In [91]:
plotTools.plotSerie(testingDataSetGP,
Title = 'GP Local Volatility Surface',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [92]:
pdeResSigmaGPTest = backtest.PDEPricerVectorized(dataSetTest, S0, nnGP, bootstrap)
pdeResSigmaGPTest.head()
Out[92]:
In [93]:
plotTools.predictionDiagnosis(pdeResSigmaGPTest,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [95]:
pdeResSigmaGPTest.to_csv(workingFolder + "pdeResSigmaGPTest.csv")
Neural Network¶
In [25]:
importlib.reload(loadData)
Out[25]:
Learning Implied volatility¶
Yon can skip the training step by loading in the left panel in colab workspace tensorflow models. These models are contained in the results folder of github repository.
In [20]:
hyperparameters = {}
#penalization coefficient
hyperparameters["lambdaLocVol"] = 0.11#1#0.001#0.01 #100
hyperparameters["lambdaSoft"] = 1#100#0.0001#10#10 #100
hyperparameters["lambdaGamma"] = 10#10000#100#10 #10000
#Derivative soft constraints parameters
hyperparameters["lowerBoundTheta"] = 0.01
hyperparameters["lowerBoundGamma"] = 0.00001
#Local variance parameters
hyperparameters["DupireVarCap"] = 10
hyperparameters["DupireVolLowerBound"] = 0.05
hyperparameters["DupireVolUpperBound"] = 0.40
#Learning scheduler coefficient
hyperparameters["LearningRateStart"] = 0.1
hyperparameters["Patience"] = 100
hyperparameters["batchSize"] = 50
hyperparameters["FinalLearningRate"] = 1e-6
hyperparameters["FixedLearningRate"] = False
#Training parameters
hyperparameters["nbUnits"] = 100 #number of units for hidden layers
hyperparameters["maxEpoch"] = 10000#10000 #maximum number of epochs
In [21]:
importlib.reload(plotTools)
Out[21]:
In [23]:
#Execute this cell if you want to fit neural network with implied volatilities
res = neuralNetwork.create_train_model_gatheral(neuralNetwork.NNArchitectureVanillaSoftGatheral,
scaledDataSet,
True,
hyperparameters,
scaler,
modelName = "convexSoftGatheralVolModel")
y_pred4G, volLocale4G, dNN_T4G, gNN_K4G, lossSerie4G = res
In [ ]:
plotTools.plotEpochLoss(lossSerie4G)
lossSerie4G.iloc[-1]
In [ ]:
y_pred4G.head()
In [43]:
# Evaluate results on the training set, you can execute that cell without training the model
res = neuralNetwork.create_eval_model_gatheral(neuralNetwork.NNArchitectureVanillaSoftGatheral,
scaledDataSet,
True,
hyperparameters,
scaler,
modelName = "convexSoftGatheralVolModel")
y_pred4G, volLocale4G, dNN_T4G, gNN_K4G, lossSerie4G = res
plotTools.modelSummaryGatheral(y_pred4G,
volLocale4G,
dNN_T4G,
gNN_K4G,
dataSet,
yMin = KMin,
yMax = KMax,
S0 = S0,
bootstrap = bootstrap,
thresholdPrice = None,
savePath = "NeuralImpliedVolTrain")
In [32]:
volLocale4G.loc[(midS0,slice(None))]
Out[32]:
In [44]:
# Evaluate results on the testing dataset, you can execute that cell without training the model
res = neuralNetwork.create_eval_model_gatheral(neuralNetwork.NNArchitectureVanillaSoftGatheral,
scaledDataSetTest,
True,
hyperparameters,
scaler,
modelName = "convexSoftGatheralVolModel")
y_pred4TestG, volLocale4TestG, dNN_T4TestG, gNN_K4TestG, lossSerie4TestG = res
plotTools.modelSummaryGatheral(y_pred4TestG,
volLocale4TestG,
dNN_T4TestG,
gNN_K4TestG,
dataSetTest,
yMin = KMin,
yMax = KMax,
S0 = S0,
bootstrap = bootstrap,
thresholdPrice = None,
savePath = "NeuralImpliedVolTest")
In [58]:
importlib.reload(plotTools)
Out[58]:
In [59]:
#Diagnosis for training results
plotTools.modelSummaryGatheral(y_pred4G,
volLocale4G,
dNN_T4G,
gNN_K4G,
dataSet,
logMoneynessScale = True,
S0 = S0,
bootstrap = bootstrap,
thresholdPrice = None,
yMin = KMin + 0.0001,
yMax = KMax)
In [60]:
#Diagnosis for testing results
plotTools.modelSummaryGatheral(y_pred4TestG,
volLocale4TestG,
dNN_T4TestG,
gNN_K4TestG,
dataSetTest,
logMoneynessScale = True,
S0 = S0,
bootstrap = bootstrap,
thresholdPrice = None,
yMin = KMin + 0.0001,
yMax = KMax)
Monte Carlo and PDE repricing backtests¶
During Monte Carlo backtest, each option in testing set is priced with an underlying which is diffused with the following SDE : $$ dS_t = \left( r_t - q_t - \frac{\sigma_{NN}^2(t, S_t)}{2} \right) dt + \sigma_{NN}(t, S_t) dW_t$$ with $\sigma_{NN}$ the neural local volatility function.
Due to computation time issue we avoid to make millions of call to neural network and we interpolate linearly neural local volatility obtained on one the two possible grid :
- the testing grid i.e. nodes $(T,K)$ of the testing set.
- an artificial grid of 10000 points to check local volatility is correctly interpolated/extrapolated. That grid is the smallest rectangle containing the minimum and maximum maturities and the minimum and maximum strikes of our dataset (union of testing and training set).
During PDE backtest, we used a crank-nicholson scheme to revaluate each option in our testing set. Time step corresponds to one day and space grid has 100 points.
In [ ]:
# Function which evaluates neural local volatility when neural network is fitted on implied volatilities
def neuralVolLocaleGatheral(s,t):
vLoc = neuralNetwork.evalVolLocaleGatheral(neuralNetwork.NNArchitectureVanillaSoftGatheral,
s, t,
dataSet,
hyperparameters,
scaler,
bootstrap,
S0,
modelName = "convexSoftGatheralVolModel")
return vLoc.dropna()
In [ ]:
#Artificial grid with 10000 points which aims to limit the impact of linear interpolation
strikeLow = min(dataSet["Strike"].min(),dataSetTest["Strike"].min()) #dataSet["Strike"].min()
strikeUp = max(dataSet["Strike"].max(),dataSetTest["Strike"].max()) #dataSet["Strike"].max()
strikeGrid = np.linspace(strikeLow, strikeUp, 100)
matLow = min(dataSet["Maturity"].min(),dataSetTest["Maturity"].min()) #dataSet["Maturity"].min()
matUp = max(dataSet["Maturity"].max(),dataSetTest["Maturity"].max()) #dataSetTest["Maturity"].max()
matGrid = np.linspace(matLow, matUp, 100)
volLocaleGrid = np.meshgrid(strikeGrid, matGrid)
In [ ]:
volLocalInterp7 = neuralVolLocaleGatheral(volLocaleGrid[0].flatten(),
volLocaleGrid[1].flatten())
volLocalInterp7.head()
In [ ]:
volLocalInterp8 = neuralVolLocaleGatheral(dataSetTest.index.get_level_values("Strike").values.flatten(),
dataSetTest.index.get_level_values("Maturity").values.flatten())
volLocalInterp8.head()
In [ ]:
#Local volatility function for backtests
nnVolLocale7 = lambda x,y : backtest.interpolatedMCLocalVolatility(volLocalInterp7, x, y)
In [ ]:
#Local volatility function for backtests
nnVolLocale8 = lambda x,y : backtest.interpolatedMCLocalVolatility(volLocalInterp8, x, y)
In [ ]:
importlib.reload(plotTools)
In [ ]:
plotTools.plotSerie(volLocalInterp7,
Title = 'Local Volatility',
az=105,
yMin=KMin,
yMax=KMax,
zAsPercent=False)
In [ ]:
plotTools.plotSerie(volLocalInterp8,
Title = 'Local Volatility',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [ ]:
nbTimeStep = 100
nbPaths = 100000
In [ ]:
#Local volatility is returned through linear interpolation on neural local volatilities obtained on the 10000 points grid.
mcResVolLocale7 = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnVolLocale7)
mcResVolLocale7.head()
In [ ]:
#Diagnose backtest results
plotTools.predictionDiagnosis(mcResVolLocale7["Price"],
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax = KMax)
In [ ]:
mcResVolLocale7.to_csv("mcResVolLocale7.csv")
In [ ]:
importlib.reload(backtest)
In [ ]:
#PDE backtest with a cranck nicholson scheme
pdeResVolLocale7 = backtest.PDEPricerVectorized(dataSetTest, S0, nnVolLocale7, bootstrap)
pdeResVolLocale7.head()
In [ ]:
#Backtest diagnosis
plotTools.predictionDiagnosis(pdeResVolLocale7,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResVolLocale7.to_csv(workingFolder + "pdeResVolLocale7.csv")
In [ ]:
In [ ]:
#Local volatility is returned through linear interpolation on neural local volatilities obtained on the testing grid.
mcResVolLocale8 = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnVolLocale8)
mcResVolLocale8.head()
In [ ]:
plotTools.predictionDiagnosis(mcResVolLocale8["Price"],
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax = KMax)
In [ ]:
mcResVolLocale8.to_csv("mcResVolLocale8.csv")
In [ ]:
pdeResVolLocale8 = backtest.PDEPricerVectorized(dataSetTest, S0, nnVolLocale8, bootstrap)
pdeResVolLocale8.head()
In [ ]:
plotTools.predictionDiagnosis(pdeResVolLocale8,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResVolLocale8.to_csv(workingFolder + "pdeResVolLocale8.csv")
NN: Learning Price¶
Yon can skip the training step by loading in the left panel in colab workspace tensorflow models. These models are contained in the results folder of github repository.
In [61]:
hyperparameters = {}
#penalization coefficient
hyperparameters["lambdaLocVol"] = 10#1#0.001#0.01 #100
hyperparameters["lambdaSoft"] = 1000#100#0.0001#10#10 #100
hyperparameters["lambdaGamma"] = 10000000#10000#100#10 #10000
#Derivative soft constraints parameters
hyperparameters["lowerBoundTheta"] = 0.01
hyperparameters["lowerBoundGamma"] = 0.00001
#Local variance parameters
hyperparameters["DupireVarCap"] = 10
hyperparameters["DupireVolLowerBound"] = 0.05
hyperparameters["DupireVolUpperBound"] = 0.40
#Learning scheduler coefficient
hyperparameters["LearningRateStart"] = 0.1
hyperparameters["Patience"] = 100
hyperparameters["batchSize"] = 50
hyperparameters["FinalLearningRate"] = 1e-6
hyperparameters["FixedLearningRate"] = False
#Training parameters
hyperparameters["nbUnits"] = 100 #number of units for hidden layers
hyperparameters["maxEpoch"] = 10000#10000 #maximum number of epochs
In [ ]:
importlib.reload(plotTools)
In [ ]:
#Execute this cell if you want to fit neural network with prices
res = neuralNetwork.create_train_model(neuralNetwork.NNArchitectureVanillaSoftDupire,
scaledDataSet,
True,
hyperparameters,
scaler,
modelName = "convexSoftVolModel")
y_pred4, volLocale4, dNN_T4, gNN_K4, lossSerie4 = res
In [ ]:
plotTools.plotEpochLoss(lossSerie4)
lossSerie4.iloc[-1]
In [ ]:
importlib.reload(plotTools)
In [62]:
# Evaluate results on the training set, you can execute that cell without training the model
res = neuralNetwork.create_eval_model(neuralNetwork.NNArchitectureVanillaSoftDupire,
scaledDataSet,
True,
hyperparameters,
scaler,
modelName = "convexSoftVolModel")
y_pred4, volLocale4, dNN_T4, gNN_K4, lossSerie4 = res
plotTools.modelSummary(y_pred4,
volLocale4,
dNN_T4,
gNN_K4,
dataSet,
S0,
bootstrap,
yMin = KMin,
yMax = KMax,
thresholdPrice = None,
removeNaN = False,
savePath = "NeuralPriceTrain")
In [ ]:
volLocale4.loc[(midS0,slice(None))]
Out[ ]:
In [63]:
# Evaluate results on the testing dataset, you can execute that cell without training the model
res = neuralNetwork.create_eval_model(neuralNetwork.NNArchitectureVanillaSoftDupire,
scaledDataSetTest,
True,
hyperparameters,
scaler,
modelName = "convexSoftVolModel")
y_pred4Test, volLocale4Test, dNN_T4Test, gNN_K4Test, lossSerie4TestG = res
plotTools.modelSummary(y_pred4Test,
volLocale4Test,
dNN_T4Test,
gNN_K4Test,
dataSetTest,
S0,
bootstrap,
yMin = KMin,
yMax = KMax,
thresholdPrice = None,
removeNaN = False,
savePath = "NeuralPriceTest")
In [ ]:
Monte Carlo and PDE repricing backtests¶
During Monte Carlo backtest, each option in testing set is priced with an underlying which is diffused with the following SDE : $$ dS_t = \left( r_t - q_t - \frac{\sigma_{NN}^2(t, S_t)}{2} \right) dt + \sigma_{NN}(t, S_t) dW_t$$ with $\sigma_{NN}$ the neural local volatility function.
Due to computation time issue we avoid to make millions of call to neural network and we interpolate linearly neural local volatility obtained on one the two possible grid :
- the testing grid i.e. nodes $(T,K)$ of the testing set.
- an artificial grid of 10000 points to check local volatility is correctly interpolated/extrapolated. That grid is the smallest rectangle containing the minimum and maximum maturities and the minimum and maximum strikes of our dataset (union of testing and training set).
During PDE backtest, we used a crank-nicholson scheme to revaluate each option in our testing set. Time step corresponds to one day and space grid has 100 points.
In [ ]:
S0
In [ ]:
# Function which evaluates neural local volatility when neural network is fitted on prices
def neuralVolLocalePrix(s,t):
vLoc = neuralNetwork.evalVolLocale(neuralNetwork.NNArchitectureVanillaSoftDupire,
s, t,
dataSet,
hyperparameters,
scaler,
bootstrap,
S0,
modelName = "convexSoftVolModel")
return vLoc.dropna()
In [ ]:
importlib.reload(neuralNetwork)
In [ ]:
#Artificial grid with 10000 points which aims to limit the impact of linear interpolation
strikeLow = min(dataSet["Strike"].min(),dataSetTest["Strike"].min()) #dataSet["Strike"].min()
strikeUp = max(dataSet["Strike"].max(),dataSetTest["Strike"].max()) #dataSet["Strike"].max()
strikeGrid = np.linspace(strikeLow, strikeUp, 100)
matLow = min(dataSet["Maturity"].min(),dataSetTest["Maturity"].min()) #dataSet["Maturity"].min()
matUp = max(dataSet["Maturity"].max(),dataSetTest["Maturity"].max()) #dataSetTest["Maturity"].max()
matGrid = np.linspace(matLow, matUp, 100)
volLocaleGrid = np.meshgrid(strikeGrid, matGrid)
In [ ]:
volLocalInterp = neuralVolLocalePrix(volLocaleGrid[0].flatten(),
volLocaleGrid[1].flatten())
volLocalInterp.head()
In [ ]:
volLocalInterp2 = neuralVolLocalePrix(dataSetTest.index.get_level_values("Strike").values.flatten(),
dataSetTest.index.get_level_values("Maturity").values.flatten())
volLocalInterp2.head()
In [ ]:
#Local volatility function for backtests
nnVolLocale = lambda x,y : backtest.interpolatedMCLocalVolatility(volLocalInterp, x, y)
In [ ]:
#Local volatility function for backtests
nnVolLocale2 = lambda x,y : backtest.interpolatedMCLocalVolatility(volLocalInterp2, x, y)
In [ ]:
plotTools.plotSerie(volLocalInterp,
Title = 'Interpolated Local Volatility Surface',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
In [ ]:
plotTools.plotSerie(volLocalInterp2,
Title = 'Interpolated Local Volatility Surface',
az=30,
yMin=KMin,
yMax=KMax,
zAsPercent=True)
Monte Carlo backtest, each option in testing set is priced with an underlying which is diffused with the following SDE : $$ dS_t = \left( r_t - q_t - \frac{\sigma_{NN}^2(t, S_t)}{2} \right) dt + \sigma_{NN}(t, S_t) dW_t$$ with $\sigma_{NN}$ the neural local volatility function.
In [ ]:
nbTimeStep = 100
nbPaths = 100000
In [ ]:
#Local volatility is returned through linear interpolation on neural local volatilities obtained on the 10000 points grid.
mcResVolLocalePrix = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnVolLocale)
mcResVolLocalePrix.head()
In [ ]:
plotTools.predictionDiagnosis(mcResVolLocalePrix["Price"],
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax = KMax)
In [ ]:
mcResVolLocalePrix.to_csv("mcResVolLocalePrix.csv")
In [ ]:
#PDE backtest with cranck-nicolson scheme
pdeResVolLocale = backtest.PDEPricerVectorized(dataSetTest, S0, nnVolLocale, bootstrap)
pdeResVolLocale.head()
In [ ]:
plotTools.predictionDiagnosis(pdeResVolLocale,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResVolLocale.to_csv(workingFolder + "pdeResVolLocale.csv")
In [ ]:
In [ ]:
In [ ]:
#Local volatility is returned through linear interpolation on neural local volatilities obtained on the testing grid.
mcResVolLocalePrix2 = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnVolLocale2)
mcResVolLocalePrix2.head()
In [ ]:
plotTools.predictionDiagnosis(mcResVolLocalePrix2["Price"],
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax = KMax)
In [ ]:
mcResVolLocalePrix2.to_csv("mcResVolLocalePrix2.csv")
In [ ]:
#PDE backtest with cranck-nicolson scheme
pdeResVolLocale2 = backtest.PDEPricerVectorized(dataSetTest, S0, nnVolLocale2, bootstrap)
pdeResVolLocale2.head()
In [ ]:
plotTools.predictionDiagnosis(pdeResVolLocale2,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResVolLocale2.to_csv(workingFolder + "pdeResVolLocale2.csv")
In [ ]:
Hyperparameter selection¶
In [ ]:
#Random selection of several hyperparameters
neuralNetwork.selectHyperparametersRandom(hyperparameters,
["lambdaLocVol","lambdaSoft","lambdaGamma"],
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
100,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = True)
In [ ]:
#hyperparameters["lambdaLocVol"] = 100
#hyperparameters["lambdaSoft"] = 100
#hyperparameters["lambdaGamma"] = 10000
hyperparameters["lambdaLocVol"] = 0.01#0.01 #100
hyperparameters["lambdaSoft"] = 0.01#10#10 #100
hyperparameters["lambdaGamma"] = 100#10 #10000
In [ ]:
#marginal selection of hypterparameters
neuralNetwork.selectHyperparameters(hyperparameters,
"lambdaLocVol",
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = True)
In [ ]:
neuralNetwork.selectHyperparameters(hyperparameters,
"DupireVarCap",
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = True)
In [ ]:
neuralNetwork.selectHyperparameters(hyperparameters,
"lambdaLocVol",
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = True)
In [ ]:
hyperparameters["lambdaLocVol"] = 100
In [ ]:
neuralNetwork.selectHyperparameters(hyperparameters,
"lambdaLocVol",
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = True)
In [ ]:
hyperparameters["nbUnits"] = 40
In [ ]:
neuralNetwork.selectHyperparameters(hyperparameters,
"nbUnits",
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = False)
In [ ]:
hyperparameters["nbUnits"] = 200
In [ ]:
neuralNetwork.selectHyperparameters(hyperparameters,
"lambdaLocVol",
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = True)
In [ ]:
hyperparameters["nbUnits"] = 40
In [ ]:
neuralNetwork.selectHyperparameters(hyperparameters,
"nbUnits",
neuralNetwork.NNArchitectureVanillaSoftGatheral,
"hyperParameters",
True,
scaledDataSet,
scaler,
trainedOnPrice = False,
logGrid = False)
In [ ]:
hyperparameters["nbUnits"] = 200
SSVI Backtesting¶
In [118]:
importlib.reload(loadData)
Out[118]:
Run SSVI Model¶
Implementation is inspired from Matlab code of Philipp Rindler :
- Philipp Rindler (2020). Gatherals and Jacquier's Arbitrage-Free SVI Volatility Surfaces (https://www.mathworks.com/matlabcentral/fileexchange/49962-gatherals-and-jacquier-s-arbitrage-free-svi-volatility-surfaces), MATLAB Central File Exchange. Retrieved November 20, 2020.
In [21]:
np.seterr(all='warn')
Out[21]:
In [22]:
importlib.reload(BS)
Out[22]:
In [23]:
#parameters, theta, maturities, pSSVI = SSVI.train_svi_surface(dataSet, S0)
SSVIModel = SSVI.SSVIModel(S0, bootstrap)
SSVIModel.fit(dataSet)
In [ ]:
In [24]:
serie = SSVIModel.eval(dataSetTest)
serieTrain = SSVIModel.eval(dataSet)
In [25]:
keptPrices = dataSet["Price"][dataSet["Price"] >= 0.0].index
plotTools.predictionDiagnosis(plotTools.selectIndex(serieTrain, keptPrices) ,
plotTools.selectIndex(dataSet["ImpliedVol"], keptPrices) ,
" Implied vol ",
yMin=KMin,
yMax=KMax)
In [26]:
impPriceTraining = plotTools.plotImpliedVolPrices(np.square(serieTrain)*serieTrain.index.get_level_values("Maturity"),
bootstrap,
S0,
dataSet,
yMin = KMin,
yMax = KMax,
thresholdPrice = None)
In [27]:
ImpVolPutTrainingSSVI = BS.vectorizedImpliedVolatilityCalibration(S0, bootstrap,
dataSet["Maturity"],
dataSet["Strike"],
dataSet["OptionType"],
impPriceTraining,
removeNaN = False)
ImpVolPutTrainingSSVI = pd.Series(ImpVolPutTrainingSSVI, index = dataSet.index).sort_index()
In [28]:
keptPrices = dataSet["Price"][dataSet["Price"] >= 0.0].index
plotTools.predictionDiagnosis(plotTools.selectIndex(ImpVolPutTrainingSSVI, keptPrices) ,
plotTools.selectIndex(dataSet["ImpliedVol"], keptPrices) ,
" Implied vol ",
yMin=KMin,
yMax=KMax)
In [29]:
keptPrices = dataSet["Price"][dataSet["Price"] >= 0.0].index
plotTools.predictionDiagnosis(plotTools.selectIndex(ImpVolPutTrainingSSVI, keptPrices),
plotTools.selectIndex(serieTrain, keptPrices),
" Implied vol ",
yMin=KMin,
yMax=KMax)
In [30]:
plotTools.predictionDiagnosis(SSVI.impliedVariance(serie),
SSVI.impliedVariance(dataSetTest["ImpliedVol"]),
" Implied total variance ",
yMin=KMin,
yMax=KMax)
In [31]:
keptPrices = dataSetTest["Price"][dataSetTest["Price"] >= 0.0].index
plotTools.predictionDiagnosis(plotTools.selectIndex(serie, keptPrices),
plotTools.selectIndex(dataSetTest["ImpliedVol"], keptPrices),
" Implied vol ",
yMin=KMin,
yMax=KMax)
In [32]:
impPriceTesting = plotTools.plotImpliedVolPrices(np.square(serie)*serie.index.get_level_values("Maturity"),
bootstrap,
S0,
dataSetTest,
yMin = KMin,
yMax = KMax,
thresholdPrice = threshold)
In [33]:
ImpVolPutTestingSSVI = BS.vectorizedImpliedVolatilityCalibration(S0, bootstrap,
dataSetTest["Maturity"],
dataSetTest["Strike"],
dataSetTest["OptionType"],
impPriceTesting,
removeNaN = False)
ImpVolPutTestingSSVI = pd.Series(ImpVolPutTestingSSVI, index = dataSetTest.index).sort_index()
In [34]:
keptPrices = dataSetTest["Price"][dataSetTest["Price"] >= 0.0].index
plotTools.predictionDiagnosis(plotTools.selectIndex(ImpVolPutTestingSSVI, keptPrices),
plotTools.selectIndex(dataSetTest["ImpliedVol"], keptPrices),
" Implied vol ",
yMin=KMin,
yMax=KMax)
In [ ]:
Estimate Local Volaitlity from SSVI¶
In [35]:
#dT, hk, dK, locVolSSVI, density = finiteDifferenceSVI(dataSet, interpolateWithSSVI)
dTTrain, hkTrain, dKTrain, locVolSSVITrain, densityTrain = SSVI.finiteDifferenceSVI(dataSet, SSVIModel.eval)
In [36]:
plotTools.diagnoseLocalVol(dTTrain,
locVolSSVITrain,
densityTrain,
SSVIModel.eval(dataSet),
dataSet,
az=320,
yMin=KMin,
yMax=KMax)
In [37]:
#dT, hk, dK, locVolSSVI, density = finiteDifferenceSVI(dataSet, interpolateWithSSVI)
dT, hk, dK, locVolSSVI, density = SSVI.finiteDifferenceSVI(dataSetTest, SSVIModel.eval)
In [38]:
plotTools.diagnoseLocalVol(dT,
locVolSSVI,
density,
SSVIModel.eval(dataSetTest),
dataSetTest,
az=320,
yMin=KMin,
yMax=KMax)
In [39]:
plotTools.saveDataModel(plotTools.removeDuplicateIndex(impPriceTraining),
plotTools.removeDuplicateIndex(locVolSSVITrain),
plotTools.removeDuplicateIndex(serieTrain) ,
"./Results/SSVITraining")
In [40]:
plotTools.saveDataModel(plotTools.removeDuplicateIndex(impPriceTesting),
plotTools.removeDuplicateIndex(locVolSSVI),
plotTools.removeDuplicateIndex(serie) ,
"./Results/SSVITesting")
In [ ]:
In [41]:
strikeLow = min(dataSet["Strike"].min(),dataSetTest["Strike"].min()) #dataSet["Strike"].min()
strikeUp = max(dataSet["Strike"].max(),dataSetTest["Strike"].max()) #dataSet["Strike"].max()
strikeGrid = np.linspace(strikeLow, strikeUp, 100)
matLow = min(dataSet["Maturity"].min(),dataSetTest["Maturity"].min()) #dataSet["Maturity"].min()
matUp = max(dataSet["Maturity"].max(),dataSetTest["Maturity"].max()) #dataSetTest["Maturity"].max()
matGrid = np.linspace(matLow, matUp, 100)
volLocaleGrid = np.meshgrid(strikeGrid, matGrid)
volLocaleGridDf = pd.DataFrame(np.vstack((np.ravel(volLocaleGrid[0]), np.ravel(volLocaleGrid[1]))).T,
columns = ["Strike", "Maturity"]).set_index(["Strike","Maturity"],drop=False)
volLocaleGridDf["ChangedStrike"] = bootstrap.changeOfVariable(volLocaleGridDf["Strike"], volLocaleGridDf["Maturity"])[0]
volLocaleGridDf["logMoneyness"] = np.log(volLocaleGridDf["ChangedStrike"] / S0)
volLocaleGridDf["OptionType"] = np.ones_like(volLocaleGridDf["logMoneyness"])
In [42]:
dT, hk, dK, locVolSSVI, density = SSVI.finiteDifferenceSVI(volLocaleGridDf, SSVIModel.eval)
In [43]:
plotTools.plotSerie(locVolSSVI.dropna(),
yMin=KMin,
yMax=KMax,
az=105,
Title = 'Local volatility')
In [44]:
plotTools.plotSerie(locVolSSVI.dropna()[locVolSSVI.dropna() <= 1.0],
yMin=KMin,
yMax=KMax,
az=105,
Title = 'Local volatility')
Backtest Local volatility SSVI¶
During Monte Carlo backtest, each option in testing set is priced with an underlying which is diffused with the following SDE : $$ dS_t = \left( r_t - q_t - \frac{\sigma_{NN}^2(t, S_t)}{2} \right) dt + \sigma_{NN}(t, S_t) dW_t$$ with $\sigma_{NN}$ the neural local volatility function.
Due to computation time issue we avoid to make millions of call to neural network and we interpolate linearly neural local volatility obtained on one the two possible grid :
- the testing grid i.e. nodes $(T,K)$ of the testing set.
- an artificial grid of 10000 points to check local volatility is correctly interpolated/extrapolated. That grid is the smallest rectangle containing the minimum and maximum maturities and the minimum and maximum strikes of our dataset (union of testing and training set).
During PDE backtest, we used a crank-nicholson scheme to revaluate each option in our testing set. Time step corresponds to one day and space grid has 100 points.
In [45]:
nbTimeStep = 100
nbPaths = 100
In [46]:
#dT, hk, dK, locVolSSVI, density = finiteDifferenceSVI(dataSet, interpolateWithSSVI)
dT, hk, dK, locVolSSVI, density = SSVI.finiteDifferenceSVI(dataSetTest, SSVIModel.eval)
nnSSVI = lambda x,y : backtest.interpolatedMCLocalVolatility(locVolSSVI, x, y)
In [47]:
strikeLow = min(dataSet["Strike"].min(),dataSetTest["Strike"].min()) #dataSet["Strike"].min()
strikeUp = max(dataSet["Strike"].max(),dataSetTest["Strike"].max()) #dataSet["Strike"].max()
strikeGrid = np.linspace(strikeLow, strikeUp, 100)
matLow = min(dataSet["Maturity"].min(),dataSetTest["Maturity"].min()) #dataSet["Maturity"].min()
matUp = max(dataSet["Maturity"].max(),dataSetTest["Maturity"].max()) #dataSetTest["Maturity"].max()
matGrid = np.linspace(matLow, matUp, 100)
volLocaleGrid = np.meshgrid(strikeGrid, matGrid)
volLocaleGridDf = pd.DataFrame(np.vstack((np.ravel(volLocaleGrid[0]), np.ravel(volLocaleGrid[1]))).T,
columns = ["Strike", "Maturity"]).set_index(["Strike","Maturity"],drop=False)
volLocaleGridDf["ChangedStrike"] = bootstrap.changeOfVariable(volLocaleGridDf["Strike"], volLocaleGridDf["Maturity"])[0]
volLocaleGridDf["logMoneyness"] = np.log(volLocaleGridDf["ChangedStrike"] / S0)
volLocaleGridDf["OptionType"] = np.ones_like(volLocaleGridDf["logMoneyness"])
In [48]:
volLocaleGridDf.head()
Out[48]:
In [ ]:
In [49]:
#dT, hk, dK, locVolSSVI, density = finiteDifferenceSVI(dataSet, interpolateWithSSVI)
dT2, hk2, dK2, locVolSSVI2, density2 = SSVI.finiteDifferenceSVI(volLocaleGridDf, SSVIModel.eval)
cleanValue = ~(dT2.isna() | density2.isna() | (dT2 < 0) | (density2 < 0) )
dT2 = dT2[cleanValue.values]
density2 = density2[cleanValue.values]
hk2 = hk2[cleanValue.values]
locVolSSVI2 = locVolSSVI2[cleanValue.values]
dK2 = dK2[cleanValue.values]
nnSSVI2 = lambda x,y : backtest.interpolatedMCLocalVolatility(locVolSSVI2, x, y)
In [ ]:
In [50]:
mcResSSVITest = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnSSVI)
mcResSSVITest.head()
Out[50]:
In [51]:
mcResSSVITest.to_csv(workingFolder + "mcResSSVITest.csv")
In [52]:
plotTools.predictionDiagnosis(mcResSSVITest["Price"],
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [53]:
mcResSSVITest2 = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnSSVI2)
mcResSSVITest2.head()
In [ ]:
plotTools.predictionDiagnosis(mcResSSVITest2["Price"],
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
mcResSSVITest2.to_csv(workingFolder + "mcResSSVITest2.csv")
In [ ]:
plotTools.plotSerie(testingDataSetSSVI[testingDataSetSSVI.index.get_level_values("Maturity") > 0.01], Title = 'SSVI Local Volatility Surface', az=30, yMin=0.0S0, yMax=2.0S0, zAsPercent=True)
In [ ]:
pdeResSigmaSSVITest = backtest.PDEPricerVectorized(dataSetTest, S0, nnSSVI, bootstrap)
pdeResSigmaSSVITest.head()
In [ ]:
plotTools.predictionDiagnosis(pdeResSigmaSSVITest,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResSigmaSSVITest.to_csv(workingFolder + "pdeResSigmaSSVITest.csv")
In [ ]:
pdeResSigmaSSVITest2 = backtest.PDEPricerVectorized(dataSetTest, S0, nnSSVI2, bootstrap)
pdeResSigmaSSVITest2.head()
In [ ]:
plotTools.predictionDiagnosis(pdeResSigmaSSVITest2,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResSigmaSSVITest2.to_csv(workingFolder + "pdeResSigmaSSVITest2.csv")
Trinomial Tree (Tikhonov)¶
MC backtest¶
During Monte Carlo backtest, each option in testing set is priced with an underlying which is diffused with the following SDE : $$ dS_t = \left( r_t - q_t - \frac{\sigma_{NN}^2(t, S_t)}{2} \right) dt + \sigma_{NN}(t, S_t) dW_t$$ with $\sigma_{NN}$ the neural local volatility function.
Due to computation time issue we avoid to make millions of call to neural network and we interpolate linearly neural local volatility obtained on one the two possible grid :
- the testing grid i.e. nodes $(T,K)$ of the testing set.
- an artificial grid of 10000 points to check local volatility is correctly interpolated/extrapolated. That grid is the smallest rectangle containing the minimum and maximum maturities and the minimum and maximum strikes of our dataset (union of testing and training set).
For the tikhonov case we only use one grid which contains the node of the trinomial tree used for the calibration algorithm. See CREPEY, Stéphane. Calibration of the local volatility in a trinomial tree using Tikhonov regularization. Inverse Problems, 2002, vol. 19, no 1, p. 91.
During PDE backtest, we used a crank-nicholson scheme to revaluate each option in our testing set. Time step corresponds to one day and space grid has 100 points.
Tikhonov calibration¶
In [ ]:
nbTimeStep = 100
nbPaths = 100000
In [122]:
nnVolLocaleTykhonov = lambda x,y : backtest.interpolatedMCLocalVolatility(localVolatility["LocalVolatility"], x, y)
In [123]:
s = testingDataSetAresky.index.get_level_values("Strike").values
t = testingDataSetAresky.index.get_level_values("Maturity").values
In [124]:
tikhonovLocVol = nnVolLocaleTykhonov(s, t)
plotTools.plotSerie(tikhonovLocVol[tikhonovLocVol.index.get_level_values("Maturity") > 0.01],
Title = 'Tikhonov Local Volatility Surface',
az=30,
yMin=0.0*S0,
yMax=2.0*S0,
zAsPercent=True)
In [125]:
mcResTikhonovTest = backtest.MonteCarloPricerVectorized(S0,
dataSetTest,
bootstrap,
nbPaths,
nbTimeStep,
nnVolLocaleTykhonov)
mcResTikhonovTest.head()
Out[125]:
In [126]:
plotTools.predictionDiagnosis(mcResTikhonovTest,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [127]:
mcResTikhonovTest.to_csv(workingFolder + "mcResTikhonovTest.csv")
In [128]:
mcResTikhonovTrain = backtest.MonteCarloPricerVectorized(S0,
dataSet,
bootstrap,
nbPaths,
nbTimeStep,
nnVolLocaleTykhonov)
mcResTikhonovTrain.head()
Out[128]:
In [129]:
plotTools.predictionDiagnosis(mcResTikhonovTrain,
dataSet["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [130]:
mcResTikhonovTrain.to_csv(workingFolder + "mcResTikhonovTrain.csv")
In [131]:
pdeResSigmaTikhonovTrain = backtest.PDEPricerVectorized(dataSet, S0, nnVolLocaleTykhonov, bootstrap)
pdeResSigmaTikhonovTrain.head()
In [ ]:
plotTools.predictionDiagnosis(pdeResSigmaTikhonovTrain,
dataSet["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResSigmaTikhonovTrain.to_csv(workingFolder + "pdeResSigmaTikhonovTrain.csv")
In [ ]:
pdeResSigmaTikhonovTest = backtest.PDEPricerVectorized(dataSetTest, S0, nnVolLocaleTykhonov, bootstrap)
pdeResSigmaTikhonovTest.head()
In [ ]:
plotTools.predictionDiagnosis(pdeResSigmaTikhonovTest,
dataSetTest["Price"],
" Price ",
yMin=KMin,
yMax=KMax)
In [ ]:
pdeResSigmaTikhonovTest.to_csv(workingFolder + "pdeResSigmaTikhonovTest.csv")
Backtest Result Comparisons¶
In [28]:
importlib.reload(backtest)
Out[28]:
In [34]:
mcResGPTest = backtest.loadMCPrices("./Results/mcResGPTest.csv",
parseHeader=0)
mcResGPTest.head()
Out[34]:
In [35]:
mcResSSVITest = backtest.loadMCPrices("./Results/mcResSSVITest.csv",
parseHeader=0)
mcResSSVITest.head()
Out[35]:
In [ ]:
backtest.rmse(mcResGPTest["Price"], mcResSSVITest["Price"])
In [39]:
pdeResSigmaGPTest = backtest.loadMCPrices("./Results/pdeResSigmaGPTest.csv")
pdeResSigmaGPTest.head()
Out[39]:
In [ ]:
Compare Implied Volatility and ignore miscalibrated implied volatilities¶
This section removes nodes where implied volatility estimation is impossible for some models (GP, NN Price). This allow a fair comparison of implied volatility for model.
In [40]:
#Detecting options which lead to implied volatilities miscalibration
gpImpVolTraining = pd.Series(pd.read_csv("./Results" + "/GPTraining.csv").set_index(["Strike","Maturity"])["ImpliedVol"],
index = plotTools.removeDuplicateIndex(dataSet).index)
gpImpVolTesting = pd.Series(pd.read_csv("./Results" + "/GPTesting.csv").set_index(["Strike","Maturity"])["ImpliedVol"],
index = plotTools.removeDuplicateIndex(dataSetTest).index)
nnImpVolTraining = pd.Series(pd.read_csv("./Results" + "/NeuralPriceTrain.csv").set_index(["Strike","Maturity"])["ImpliedVol"],
index = plotTools.removeDuplicateIndex(dataSet).index)
nnImpVolTesting = pd.Series(pd.read_csv("./Results" + "/NeuralPriceTest.csv").set_index(["Strike","Maturity"])["ImpliedVol"],
index = plotTools.removeDuplicateIndex(dataSetTest).index)
In [41]:
gpImpVolTrainingError = gpImpVolTraining[gpImpVolTraining <= 1e-9].append(gpImpVolTraining[gpImpVolTraining >= 2])
nnImpVolTrainingError = nnImpVolTraining[nnImpVolTraining <= 1e-9].append(nnImpVolTraining[nnImpVolTraining >= 2])
trainingFilterImpVol = gpImpVolTrainingError.append(nnImpVolTrainingError).index.unique()
In [42]:
gpImpVolTestingError = gpImpVolTesting[gpImpVolTesting <= 1e-9].append(gpImpVolTesting[gpImpVolTesting >= 2])
nnImpVolTestingError = nnImpVolTesting[nnImpVolTesting <= 1e-9].append(nnImpVolTesting[nnImpVolTesting >= 2])
testingFilterImpVol = gpImpVolTestingError.append(nnImpVolTestingError).index.unique()
In [ ]:
In [43]:
#Load results from another model
dfModelTraining = pd.read_csv("./Results" + "/NeuralImpliedVolTrain.csv").set_index(["Strike","Maturity"])
dfModelTraining = pd.DataFrame(dfModelTraining.values,
index = plotTools.removeDuplicateIndex(dataSet).index,
columns = dfModelTraining.columns)
In [44]:
keptVols = dataSet["ImpliedVol"].drop(trainingFilterImpVol).index
plotTools.predictionDiagnosis(plotTools.selectIndex(dfModelTraining["ImpliedVol"], keptVols) ,
plotTools.selectIndex(dataSet["ImpliedVol"], keptVols) ,
" Implied vol ",
yMin=KMin,
yMax=KMax)
In [45]:
dfModelTesting = pd.read_csv("./Results" + "/NeuralImpliedVolTest.csv").set_index(["Strike","Maturity"])
dfModelTesting = pd.DataFrame(dfModelTesting.values,
index = plotTools.removeDuplicateIndex(dataSetTest).index,
columns = dfModelTesting.columns)
In [46]:
keptVols = dataSetTest["ImpliedVol"].drop(testingFilterImpVol).index
plotTools.predictionDiagnosis(plotTools.selectIndex(dfModelTesting["ImpliedVol"], keptVols) ,
plotTools.selectIndex(dataSetTest["ImpliedVol"], keptVols) ,
" Implied vol ",
yMin=KMin,
yMax=KMax)
